
# Bernd Beber and Chris Blattman, "The Logic of Child Soldiering and Coercion"
# Equilibrium solution
# Last modified February 2012
# R version 2.14.1

# Load package needed for figures
    library(splines)

# Leader's utility
    uL <- function(theta, a, uG, uB, k) {
        theta * a * (1 - uG^2) + (1 - a * theta) * (-ifelse(uB<0, k, 1) * uB^2)
        }

# Recruit's utility
    uR <- function(theta, a, uG, uB, n) {
        theta * a * uG + (1 - a * theta) * uB - theta^n * a^2
        }

# Function to locate equilibria
    eq <- function(y, n, k, ul, uh, a.min=0, a.max=1) {
        # Set up list for valid candidate solutions
        solutions <- c()
        # Participation constraint binds
        uB <- function(u, theta, a, n=n) u - theta^n * a^2
        uG <- function(u, theta, a, n=n) u - theta^n * a^2 + 2 * theta^(n-1) * a
        lambda1 <- function(u, theta, k, a, n=n) 2 * a * (1 - a * theta) * ((u - theta^n * a^2) * (1 - k) + 2 * theta^(n-1) * a)
        lambda2 <- function(u, theta, k, a, n=n) 2 * (u - theta^n * a^2) * ((1 - a * theta) * k + a * theta) + 4 * theta^n * a^2
        max.a <- function(a, u, theta, k, n=n) {
            .uB <- uB(u=u, theta=theta, a=a, n=n)
            .uG <- uG(u=u, theta=theta, a=a, n=n)
            .lambda1 <- lambda1(u=u, theta=theta, k=k, a=a, n=n)
            .lambda2 <- lambda2(u=u, theta=theta, k=k, a=a, n=n)
            theta * a * (1 - .uG^2) + (1 - theta * a) * (-k * .uB^2) - .lambda1 * (2 * theta^n * a - theta * (.uG - .uB)) - .lambda2 * (u - theta * a * .uG - (1 - theta * a) * .uB + theta^n * a^2)
        }
        op <- optimize(max.a, c(a.min, a.max), u=ul, theta=y, k=k, n=n, maximum=T, tol=.Machine$double.eps)$maximum
        .alneg <- ifelse(uB(u=ul, theta=y, a=op, n=n) <= 0 &
                         uL(theta=y, a=op, uG=uG(u=ul, theta=y, a=op, n=n), uB=uB(u=ul, theta=y, a=op, n=n), k=k) >= 0 &
                         lambda1(u=ul, theta=y, k=k, a=op, n=n) >= 0 &
                         lambda2(u=ul, theta=y, k=k, a=op, n=n) >= 0, op, NA)
        if(!is.na(.alneg)) solutions <- rbind(solutions,
                                              data.frame(theta=y, a=.alneg, uG=uG(u=ul, theta=y, a=.alneg, n=n), uB=uB(u=ul, theta=y, a=.alneg, n=n),
                                                         k=k, n=n, ul=ul, uh=uh, y=y,
                                                         uL=uL(theta=y, a=.alneg, uG=uG(u=ul, theta=y, a=.alneg, n=n), uB=uB(u=ul, theta=y, a=.alneg, n=n), k=k),                                              
                                                         uR=uR(theta=y, a=.alneg, uG=uG(u=ul, theta=y, a=.alneg, n=n), uB=uB(u=ul, theta=y, a=.alneg, n=n), n=n),
                                                         lambda1=lambda1(u=ul, theta=y, k=k, a=.alneg, n=n), lambda2=lambda2(u=ul, theta=y, k=k, a=.alneg, n=n),
                                                         case="PC binds", effort=ifelse(isTRUE(all.equal(.alneg, 1, tol=.001)), "Upper bound", "Interior"), type="Child"))
        op <- optimize(max.a, c(a.min, a.max), u=ul, theta=y, k=1, n=n, maximum=T, tol=.Machine$double.eps)$maximum
        .alpos <- ifelse(uB(u=ul, theta=y, a=op, n=n) >= 0 &
                         uL(theta=y, a=op, uG=uG(u=ul, theta=y, a=op, n=n), uB=uB(u=ul, theta=y, a=op, n=n), k=k) >= 0 &
                         lambda1(u=ul, theta=y, k=1, a=op, n=n) >= 0 &
                         lambda2(u=ul, theta=y, k=1, a=op, n=n) >= 0, op, NA)
        if(!is.na(.alpos)) solutions <- rbind(solutions,
                                              data.frame(theta=y, a=.alpos, uG=uG(u=ul, theta=y, a=.alpos, n=n), uB=uB(u=ul, theta=y, a=.alpos, n=n),
                                                         k=k, n=n, ul=ul, uh=uh, y=y,
                                                         uL=uL(theta=y, a=.alpos, uG=uG(u=ul, theta=y, a=.alpos, n=n), uB=uB(u=ul, theta=y, a=.alpos, n=n), k=k),                                              
                                                         uR=uR(theta=y, a=.alpos, uG=uG(u=ul, theta=y, a=.alpos, n=n), uB=uB(u=ul, theta=y, a=.alpos, n=n), n=n),
                                                         lambda1=lambda1(u=ul, theta=y, k=k, a=.alpos, n=n), lambda2=lambda2(u=ul, theta=y, k=k, a=.alpos, n=n),
                                                         case="PC binds", effort=ifelse(isTRUE(all.equal(.alpos, 1, tol=.001)), "Upper bound", "Interior"), type="Child"))
        op <- optimize(max.a, c(a.min, a.max), u=uh, theta=1, k=k, n=n, maximum=T, tol=.Machine$double.eps)$maximum
        .ahneg <- ifelse(uB(u=uh, theta=1, a=op, n=n) <= 0 &
                         uL(theta=1, a=op, uG=uG(u=uh, theta=1, a=op, n=n), uB=uB(u=uh, theta=1, a=op, n=n), k=k) >= 0 &
                         lambda1(u=uh, theta=1, k=k, a=op, n=n) >= 0 &
                         lambda2(u=uh, theta=1, k=k, a=op, n=n) >= 0, op, NA)
        if(!is.na(.ahneg)) solutions <- rbind(solutions,
                                              data.frame(theta=1, a=.ahneg, uG=uG(u=uh, theta=1, a=.ahneg, n=n), uB=uB(u=uh, theta=1, a=.ahneg, n=n),
                                                         k=k, n=n, ul=ul, uh=uh, y=y,
                                                         uL=uL(theta=1, a=.ahneg, uG=uG(u=uh, theta=1, a=.ahneg, n=n), uB=uB(u=uh, theta=1, a=.ahneg, n=n), k=k),                                              
                                                         uR=uR(theta=1, a=.ahneg, uG=uG(u=uh, theta=1, a=.ahneg, n=n), uB=uB(u=uh, theta=1, a=.ahneg, n=n), n=n),
                                                         lambda1=lambda1(u=uh, theta=1, k=k, a=.ahneg, n=n), lambda2=lambda2(u=uh, theta=1, k=k, a=.ahneg, n=n),
                                                         case="PC binds", effort=ifelse(isTRUE(all.equal(.ahneg, 1, tol=.001)), "Upper bound", "Interior"), type="Adult"))
        op <- optimize(max.a, c(a.min, a.max), u=uh, theta=1, k=1, n=n, maximum=T, tol=.Machine$double.eps)$maximum
        .ahpos <- ifelse(uB(u=uh, theta=1, a=op, n=n) >= 0 &
                         uL(theta=1, a=op, uG=uG(u=uh, theta=1, a=op, n=n), uB=uB(u=uh, theta=1, a=op, n=n), k=k) >= 0 &
                         lambda1(u=uh, theta=1, k=1, a=op, n=n) >= 0 &
                         lambda2(u=uh, theta=1, k=1, a=op, n=n) >= 0, op, NA)
        if(!is.na(.ahpos)) solutions <- rbind(solutions,
                                              data.frame(theta=1, a=.ahpos, uG=uG(u=uh, theta=1, a=.ahpos, n=n), uB=uB(u=uh, theta=1, a=.ahpos, n=n),
                                                         k=k, n=n, ul=ul, uh=uh, y=y,
                                                         uL=uL(theta=1, a=.ahpos, uG=uG(u=uh, theta=1, a=.ahpos, n=n), uB=uB(u=uh, theta=1, a=.ahpos, n=n), k=k),                                              
                                                         uR=uR(theta=1, a=.ahpos, uG=uG(u=uh, theta=1, a=.ahpos, n=n), uB=uB(u=uh, theta=1, a=.ahpos, n=n), n=n),
                                                         lambda1=lambda1(u=uh, theta=1, k=k, a=.ahpos, n=n), lambda2=lambda2(u=uh, theta=1, k=k, a=.ahpos, n=n),
                                                         case="PC binds", effort=ifelse(isTRUE(all.equal(.ahpos, 1, tol=.001)), "Upper bound", "Interior"), type="Adult"))
        # Participation constraint doesn't bind
        uB <- function(theta, a, n=n, k=k) (-4 * theta^n * a^2) / (theta + 2 * k * (1 - a * theta))
        uG <- function(theta, a, n=n, k=k) 2 * theta^(n-1) * a + (-4 * theta^n * a^2) / (theta + 2 * k * (1 - a * theta))
        lambda1 <- function(theta, a, n=n, k=k) 2 * a * uG(theta=theta, a=a, n=n, k=k)
        max.a <- function(a, theta, k, n=n) {
            .uB <- uB(theta=theta, a=a, n=n, k=k)
            .uG <- uG(theta=theta, a=a, n=n, k=k)
            .lambda1 <- lambda1(theta=theta, a=a, n=n, k=k)
            theta * a * (1 - .uG^2) + (1 - theta * a) * (-k * .uB^2) - .lambda1 * (2 * theta^n * a - theta * (.uG - .uB))
        }
        op <- optimize(max.a, c(a.min, a.max), theta=y, k=k, n=n, maximum=T, tol=.Machine$double.eps)$maximum
        .al <- ifelse(uL(theta=y, a=op, uG=uG(theta=y, a=op, n=n, k=k), uB=uB(theta=y, a=op, n=n, k=k), k=k) >= 0 &
                      lambda1(theta=y, a=op, n=n, k=k) >= 0 &
                      uR(theta=y, a=op, uG=uG(theta=y, a=op, n=n, k=k), uB=uB(theta=y, a=op, n=n, k=k), n=n) >= ul, op, NA)
        if(!is.na(.al)) solutions <- rbind(solutions,
                                           data.frame(theta=y, a=.al, uG=uG(theta=y, a=.al, n=n, k=k), uB=uB(theta=y, a=.al, n=n, k=k),
                                                      k=k, n=n, ul=ul, uh=uh, y=y,
                                                      uL=uL(theta=y, a=.al, uG=uG(theta=y, a=.al, n=n, k=k), uB=uB(theta=y, a=.al, n=n, k=k), k=k),                                              
                                                      uR=uR(theta=y, a=.al, uG=uG(theta=y, a=.al, n=n, k=k), uB=uB(theta=y, a=.al, n=n, k=k), n=n),
                                                      lambda1=lambda1(theta=y, a=.al, n=n, k=k), lambda2=0,
                                                      case="PC does not bind", effort=ifelse(isTRUE(all.equal(.al, 1, tol=.001)), "Upper bound", "Interior"), type="Child"))
        op <- optimize(max.a, c(a.min, a.max), theta=1, k=k, n=n, maximum=T, tol=.Machine$double.eps)$maximum
        .ah <- ifelse(uL(theta=1, a=op, uG=uG(theta=1, a=op, n=n, k=k), uB=uB(theta=1, a=op, n=n, k=k), k=k) >= 0 &
                      lambda1(theta=1, a=op, n=n, k=k) >= 0 &
                      uR(theta=1, a=op, uG=uG(theta=1, a=op, n=n, k=k), uB=uB(theta=1, a=op, n=n, k=k), n=n) >= uh, op, NA)
        if(!is.na(.ah)) solutions <- rbind(solutions,
                                           data.frame(theta=1, a=.ah, uG=uG(theta=1, a=.ah, n=n, k=k), uB=uB(theta=1, a=.ah, n=n, k=k),
                                                      k=k, n=n, ul=ul, uh=uh, y=y,
                                                      uL=uL(theta=1, a=.ah, uG=uG(theta=1, a=.ah, n=n, k=k), uB=uB(theta=1, a=.ah, n=n, k=k), k=k),                                              
                                                      uR=uR(theta=1, a=.ah, uG=uG(theta=1, a=.ah, n=n, k=k), uB=uB(theta=1, a=.ah, n=n, k=k), n=n),
                                                      lambda1=lambda1(theta=1, a=.ah, n=n, k=k), lambda2=0,
                                                      case="PC does not bind", effort=ifelse(isTRUE(all.equal(.ah, 1, tol=.001)), "Upper bound", "Interior"), type="Adult"))
        # Corner solutions, no effort
        # Need non-positive reservation value to meet participation constraint
        if(uh <= 0) solutions <- rbind(solutions,
                                       data.frame(theta=1, a=0, uG=0, uB=0,  # Can set uG to anything, since P(G)=0
                                                  k=k, n=n, ul=ul, uh=uh, y=y,
                                                  uL=0, uR=0, lambda1=0, lambda2=0,
                                                  case="No production", effort="Lower bound", type="Adult"))
        if(ul <= 0) solutions <- rbind(solutions,
                                       data.frame(theta=y, a=0, uG=0, uB=0,  # Can set uG to anything, since P(G)=0
                                                  k=k, n=n, ul=ul, uh=uh, y=y,
                                                  uL=0, uR=0, lambda1=0, lambda2=0,
                                                  case="No production", effort="Lower bound", type="Child"))
        # Order solutions by leader's utility
        if(length(solutions)>0) solutions <- solutions[order(solutions[,"uL"], decreasing=T), ]
        return(solutions)
    }

# Function to return smoothed values from function �f� for input values �pred�
# where smoother is modeled on function values over input �grid� using a spline with degrees of freedom �df�
# Requires package �splines�
    sm <- function(f, pred, grid, df=15) {
        .f <- Vectorize(f)
        .f.out <- cbind(x=grid, y=.f(grid))
        # Alternative for smoothing:
        # fit <- lm(y ~ poly(x, 5), as.data.frame(.f.out))
        fit <- lm(y ~ ns(x, df), as.data.frame(.f.out))
        predict(fit, data.frame(x=pred))
    }

# FIGURE 1:
# Leader's utility, by recruit type and punishment cost
    postscript("Figure1.eps", width=7, height=7, bg="white", paper="special", horizontal=F, onefile=F)
        # Set up plot
        plot(NA, xlim=c(.01, 1), ylim=c(0, .4), ylab="Leader's utility", xlab="Productivity parameter y")
        # Utility from child recruitment, k=1
        uL.y <- function(x) {
            tmp <- eq(y=x, n=0, k=1, ul=-.2, uh=-.2)
            tmp <- tmp[tmp[,"type"]=="Child", "uL"]
            ifelse(length(tmp)>0, max(tmp), NA) }
        ys <- sm(uL.y, seq(.01, 1, l=100), seq(.01, 1, l=1000))
        lines(seq(.01, 1, l=100), c(uL.y(.01), ys[-c(1,100)], uL.y(1)))
        # Utility from adult recruitment, k=1
        uL.y <- function(x) {
            tmp <- eq(y=x, n=0, k=1, ul=-.2, uh=-.2)
            tmp <- tmp[tmp[,"type"]=="Adult", "uL"]
            ifelse(length(tmp)>0, max(tmp), NA) }
        V.uL.y <- Vectorize(uL.y)
        curve(V.uL.y, add=T)
        # Utility from child recruitment, k=.05
        uL.y <- function(x) {
            tmp <- eq(y=x, n=0, k=.05, ul=-.2, uh=-.2)
            tmp <- tmp[tmp[,"type"]=="Child", "uL"]
            ifelse(length(tmp)>0, max(tmp), NA) }
        ys <- sm(uL.y, seq(.01, 1, l=100), seq(.01, 1, l=1000))
        lines(seq(.01, 1, l=100), c(uL.y(.01), ys[-c(1,100)], uL.y(1)))
        # Utility from adult recruitment, k=.05
        uL.y <- function(x) {
            tmp <- eq(y=x, n=0, k=.05, ul=-.2, uh=-.2)
            tmp <- tmp[tmp[,"type"]=="Adult", "uL"]
            ifelse(length(tmp)>0, max(tmp), NA) }
        V.uL.y <- Vectorize(uL.y)
        curve(V.uL.y, add=T)
        # Add labels
        text(.6, .048, "Child soldiering", pos=1, cex=.9)
        text(.6, .275, "Adult recruitment", pos=1, cex=.9)
        text(.6, .375, "Adult recruitment", pos=1, cex=.9)
        text(0, .338, "k = 1", pos=4, cex=.9)
        text(0, .238, "k = .05", pos=4, cex=.9)
    dev.off()

# FIGURE 2:
# Leader's utility when indoctrination is effective
    postscript("Figure2.eps", width=7, height=7, bg="white", paper="special", horizontal=F, onefile=F)
        plot(NA, xlim=c(0, 8), ylim=c(0, .8), ylab="Leader's utility", xlab="Indoctrination parameter n")
        # Utility from child recruitment, y=.7
        uL.n <- function(x) {
            tmp <- eq(y=.7, n=x, k=1, ul=0, uh=0)
            tmp <- tmp[tmp[,"type"]=="Child", "uL"]
            ifelse(length(tmp)>0, max(tmp), NA) }
        ys <- sm(uL.n, seq(0, 8, l=100), seq(0, 8, l=1000))
        lines(seq(0, 8, l=100), ys)
        # Utility from child recruitment, y=.4
        uL.n <- function(x) {
            tmp <- eq(y=.4, n=x, k=1, ul=0, uh=0)
            tmp <- tmp[tmp[,"type"]=="Child", "uL"]
            ifelse(length(tmp)>0, max(tmp), NA) }
        ys <- sm(uL.n, seq(0, 8, l=100), seq(0, 8, l=1000))
        lines(seq(0, 8, l=100), ys)
        # Utility from adult recruitment, invariant in y
        uL.n <- function(x) {
            tmp <- eq(y=.7, n=x, k=1, ul=0, uh=0)
            tmp <- tmp[tmp[,"type"]=="Adult", "uL"]
            ifelse(length(tmp)>0, max(tmp), NA) }
        V.uL.n <- Vectorize(uL.n)
        curve(V.uL.n, add=T)
        # Add labels
        text(8, .625, "Child soldiering, y = .7", pos=2, cex=.9)
        text(8, .37, "Child soldiering, y = .4", pos=2, cex=.9)
        text(8, .195, "Adult recruitment", pos=2, cex=.9)
    dev.off()

# FIGURE 3:
# The effectiveness of indoctrination as punishment becomes cheap
    postscript("Figure3.eps", width=7, height=7, bg="white", paper="special", horizontal=F, onefile=F)
        par(mfrow=c(2,2))
        par(mar=c(5, 4, 2, 2) + .1)
        uL.n.k <- function(k, label.x=F, label.y=F, label.lines=F) {
            plot(NA, xlim=c(0,2), ylim=c(0,.4), ylab="", xlab="", main=paste("k =", k), cex.main=1)
            # Utility from child recruitment, y=.9
            uL.n <- function(x) {
                tmp <- eq(y=.9, n=x, k=k, ul=-.2, uh=0)
                tmp <- tmp[tmp[,"type"]=="Child", "uL"]
                ifelse(length(tmp)>0 & max(tmp)>0, max(tmp), NA) }
            ys <- sm(uL.n, seq(0, 1.91, l=100), seq(0, 1.91, l=1000))
            lines(seq(0, 1.91, l=100), ys)
            # Utility from child recruitment, y=.7
            uL.n <- function(x) {
                tmp <- eq(y=.7, n=x, k=k, ul=-.2, uh=0)
                tmp <- tmp[tmp[,"type"]=="Child", "uL"]
                ifelse(length(tmp)>0 & max(tmp)>0, max(tmp), NA) }
            ys <- sm(uL.n, seq(0, 1.91, l=100), seq(0, 1.91, l=1000))
            lines(seq(0, 1.91, l=100), ys)
            # Utility from child recruitment, y=.5
            uL.n <- function(x) {
                tmp <- eq(y=.5, n=x, k=k, ul=-.2, uh=0)
                tmp <- tmp[tmp[,"type"]=="Child", "uL"]
                ifelse(length(tmp)>0 & max(tmp)>0, max(tmp), NA) }
            ys <- sm(uL.n, seq(0, 1.91, l=100), seq(0, 1.91, l=1000))
            lines(seq(0, 1.91, l=100), ys)
            # Utility from adult recruitment, invariant in y
            uL.n <- function(x) {
                tmp <- eq(y=.7, n=x, k=k, ul=-.2, uh=0)
                tmp <- tmp[tmp[,"type"]=="Adult", "uL"]
                ifelse(length(tmp)>0 & max(tmp)>0, max(tmp), NA) }
            V.uL.n <- Vectorize(uL.n)
            curve(V.uL.n, xlim=c(0,1.91), add=T)
            # Add labels
            if(label.x) text(-.53, -.14, "Indoctrination parameter n", xpd=NA)
            if(label.y) text(-.55, -.123, "Leader's utility", xpd=NA, srt=90)
            if(label.lines) {
                text(.6, .207, "y = .9", pos=1, cex=.9)
                text(.6, .135, "y = .7", pos=1, cex=.9)
                text(.6, .081, "y = .5", pos=1, cex=.9)
            }
        }
        # Repeat for different values of k
        uL.n.k(1, label.y=T, label.lines=T)
        uL.n.k(.5)
        uL.n.k(.2)
        uL.n.k(.05, label.x=T)
    dev.off()

# FIGURE 4:
# Leader's utility when adult-child reservation values differ
    postscript("Figure4.eps", width=7, height=7, bg="white", paper="special", horizontal=F, onefile=F)
        plot(NA, xlim=c(0, .8), ylim=c(0, .2), ylab="Leader's utility", xlab=expression(paste("Distance between reservation values, with constant ", underline(u)[H])))
        # Utility from child recruitment, y=.9
        uL.delta <- function(x) {
            tmp <- eq(y=.9, n=0, k=1, ul=.3-x, uh=.3)
            tmp <- tmp[tmp[,"type"]=="Child", "uL"]
            ifelse(length(tmp)>0, max(tmp), NA) }
        ys <- sm(uL.delta, seq(0, .8, l=100), c(seq(0, .39, l=500), seq(.41, .8, l=500)))
        lines(seq(0, .8, l=100), ys)
        # Utility from child recruitment, y=.7
        uL.delta <- function(x) {
            tmp <- eq(y=.7, n=0, k=1, ul=.3-x, uh=.3)
            tmp <- tmp[tmp[,"type"]=="Child", "uL"]
            ifelse(length(tmp)>0, max(tmp), NA) }
        ys <- sm(uL.delta, seq(0, .8, l=100), seq(0, .8, l=1000))
        lines(seq(0, .8, l=100)[ys>0], ys[ys>0])
        # Utility from child recruitment, y=.4
        uL.delta <- function(x) {
            tmp <- eq(y=.4, n=0, k=1, ul=.3-x, uh=.3)
            tmp <- tmp[tmp[,"type"]=="Child", "uL"]
            ifelse(length(tmp)>0, max(tmp), NA) }
        ys <- sm(uL.delta, seq(.1, .8, l=100), seq(.1, .8, l=1000))
        lines(seq(.1, .8, l=100)[ys>0], ys[ys>0])
        # Utility from adult recruitment, invariant to y
        uL.delta <- function(x) {
            tmp <- eq(y=.7, n=0, k=1, ul=.3-x, uh=.3)
            tmp <- tmp[tmp[,"type"]=="Adult", "uL"]
            ifelse(length(tmp)>0, max(tmp), NA) }
        V.uL.delta <- Vectorize(uL.delta)
        curve(V.uL.delta, add=T)
        # Add labels
        text(.8, .189, "Child soldiering, y = .9", pos=2, cex=.9)
        text(.8, .108, "Child soldiering, y = .7", pos=2, cex=.9)
        text(.8, .078, "Adult recruitment", pos=2, cex=.9)
        text(.8, .037, "Child soldiering, y = .4", pos=2, cex=.9)
    dev.off()

# FIGURE 5:
# The impact of differing reservation utilities as punishment becomes cheap
    postscript("Figure5.eps", width=7, height=7, bg="white", paper="special", horizontal=F, onefile=F)
        par(mfrow=c(2,2))
        par(mar=c(5, 4, 2, 2) + .1)
        uL.delta.k <- function(k, label.x=F, label.y=F, label.lines=F) {
            plot(NA, xlim=c(0, .8), ylim=c(0, .2), ylab="", xlab="", main=paste("k =", k), cex.main=1)
            # Utility from child recruitment, y=.9
            uL.delta <- function(x) {
                tmp <- eq(y=.9, n=0, k=k, ul=.3-x, uh=.3)
                tmp <- tmp[tmp[,"type"]=="Child", "uL"]
                ifelse(length(tmp)>0, max(tmp), NA) }
            ys <- sm(uL.delta, seq(0, .8, l=100), c(seq(0, .39, l=500), seq(.41, .8, l=500)))
            lines(seq(0, .8, l=100), ys)
            # Utility from child recruitment, y=.7
            uL.delta <- function(x) {
                tmp <- eq(y=.7, n=0, k=k, ul=.3-x, uh=.3)
                tmp <- tmp[tmp[,"type"]=="Child", "uL"]
                ifelse(length(tmp)>0, max(tmp), NA) }
            ys <- sm(uL.delta, seq(0, .8, l=100), c(seq(0, .41, l=500), seq(.48, .8, l=500)))
            lines(seq(0, .8, l=100)[ys>0], ys[ys>0])
            # Utility from child recruitment, y=.5
            uL.delta <- function(x) {
                tmp <- eq(y=.5, n=0, k=k, ul=.3-x, uh=.3)
                tmp <- tmp[tmp[,"type"]=="Child", "uL"]
                ifelse(length(tmp)>0, max(tmp), NA) }
            ys <- sm(uL.delta, seq(.05, .8, l=100), c(seq(.05, .4, l=500), seq(.55, .8, l=500)))
            lines(seq(.05, .8, l=100)[ys>0], ys[ys>0])
            # Utility from adult recruitment, invariant in y
            uL.delta <- function(x) {
                tmp <- eq(y=.7, n=0, k=k, ul=.3-x, uh=.3)
                tmp <- tmp[tmp[,"type"]=="Adult", "uL"]
                ifelse(length(tmp)>0, max(tmp), NA) }
            V.uL.delta <- Vectorize(uL.delta)
            curve(V.uL.delta, add=T)
            # Add labels
            if(label.x) text(-.21, -.07, expression(paste("Distance between reservation values, with constant ", underline(u)[H])), xpd=NA)
            if(label.y) text(-.22, -.0615, "Leader's utility", xpd=NA, srt=90)
            if(label.lines) {
                text(.8, .174, "y = .9", pos=2, cex=.9)
                text(.8, .093, "y = .7", pos=2, cex=.9)
                text(.8, .04, "y = .5", pos=2, cex=.9)
            }
            # Clean up frame
            polygon(c(0, .8, .8, 0), c(.2, .2, .3, .3), border=NA, col="white")
            box()
        }
        uL.delta.k(1, label.y=T, label.lines=T)
        uL.delta.k(.5)
        uL.delta.k(.2)
        uL.delta.k(.05, label.x=T)
    dev.off()
